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investigate  the  sputtering  yields  from  both  liquid  and  solid 
metal  targets.  The  system  consisted  of  1.0  kev  Argon  ions 
bombarding  Rhodium  targets.  The  embedded  atom  method  of 
calculating  potentials  was  used  with  a  modified 
Moliere/Morse  potential  function.  The  yields  from  the  solid 
and  liquid  targets  were  compared  with  the  liquid  showing  a 
slightly  higher  yield  than  the  solid.  The  liquid  was 
simulated  by  random  displacements  of  the  atoms  from  a  solid 
crystal  lattice.  Changing  the  seed,  used  by  the  random 
number  generator  to  produce  the  liquid,  effected  the 
sputtering  yield  similar  to  moving  the  impact  point.  Four 
different  sampling  methods  were  investigated  which  produced 
similar  results.  Hence,  the  models  described  in  this  thesis 
should  provide  a  basis  for  general  sputtering  simulations  of 
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I .  INTRODUCTION 


A.  BACKGROUND  ON  SPUTTERING 

Sputtering  is  the  ejection  of  atoms  or  molecules  from 
the  surface  of  a  target  material  as  a  result  of  the 
bombardment  of  that  surface  by  particles.  The  general  area 
of  surface  bombardment  is  usually  divided  into  the  two 
categories,  of  forward  and  backward  effects.  The  formation 
of  vacancies  and  interstitial  atoms  is  an  example  of  a 
forward  effect  which  can  occur;  a  vacancy  and  an 
interstitial  atom  together  making  up  a  Frenkel  Pair. 
Scattering  of  the  incident  particles,  ejection  of  electrons, 
and  ejection  of  photons  are  just  a  few  examples  of  the  many 
backward  processes  that  can  occur. 

The  bombarding  particles  can  be  neutral  or  ionized  atoms 
or  molecules,  neutrons,  electrons,  or  photons.  The  ejected 
target  particles  can  be  neutral  or  ionized,  in  much  of  the 
research  that  has  been  done,  the  incident  particles  have 
been  ions,  and  the  target  materials  studied  have  been  atoms, 
so  that  unless  a  specific  differentiation  is  needed,  the 
terms  'incident  ions'  and  'ejected  atoms'  will  be  commonly 
used. 

In  ion  bombardment  there  are  a  number  of  important  areas 
that  require  detailed  study.  In  ion  implantation  the  state 
of  the  impinging  ion  needs  to  be  determined.  How  does  it 
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lose  its  energy,  what  is  its  path,  where  does  it  come  to 

rest?  In  dynamic  recoil  mixing,  it  is  the  state  of  the 

first  atom  hit  by  the  impinging  ion.  This  atom  is  called 

the  primary  knock-on  atom  ( PKA) .  In  the  case  where  more 

than  one  atom  is  hit  by  the  impinging  ion,  multiple  PKA's  ,  ] 

can  be  produced.  In  sputtering,  all  the  atoms  of  the  target  ^ 

that  are  involved  in  the  interaction,  can  potentially  be 

ejected.  These  moving  atoms  form  the  collision  cascade. 

The  bombardment  processes,  both  forward  or  backward,  can  be 
described  by  the  motion  of  these  atoms. 

The  quantitative  measure  of  surface  erosion  by 
sputtering  is  the  yield.  The  yield  is  defined  as  the  number 
of  ejected  atoms  per  incident  ion.  Researchers  have 
investigated  the  dependence  of  the  yield  on  physical 
quantities  such  as  the  incident  ion’s  energy,  species, 
charge,  and  angle  of  incidence;  and  the  target's  composition 
and  structure. 

Sputtering  can  be  categorized  according  to  the 
microscopic  mechanisms  that  cause  the  ejection  of  atoms  from 
the  target  material.  These  include  collisional,  thermal, 
chemical,  and  electronic  processes.  Only  collisional 
sputtering  will  be  investigated  in  this  project. 

A  bombardment  begins  as  a  collision  between  an  incident 
ion  and  a  target  atom.  If  enough  energy  is  transferred  to 
the  target  atom,  enabling  it  to  overcome  the  binding  energy 
of  its  lattice  position,  a  recoil  atom  is  created  which  can 
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collide  with  other  target  atoms  causing  a  cascade.  When  the 
energy  transferred  back  to  a  surface  atom  is  greater  than 
the  binding  energy  at  the  surface,  the  atom  is  ejected. 
This  energy  is  called  the  sputtering  threshold. 

The  interaction  that  occurs  between  the  ion  and  atom  or 
between  the  recoil  atom  and  subsequent  target  atoms  can  be 
described  using  the  classical  mechanics  of  two-body  elastic 
collisions.  If  electronic  energy  loss  is  ignored,  energy 
and  momentum  are  conserved.  We  look  specifically  at  the 
transfer  of  energy  and  momentum  from  the  impinging  ion  to 
the  atoms  in  the  target.  Inelastic  energy  losses  such  as 
electron  excitation  can  occur  as  part  of  the  dynamics,  but 
consideration  of  these  processes  is  beyond  the  scope  of  this 
project.  In  order  to  carry  out  the  calculations  it  is 
necessary  to  specify  the  potential  and  the  forces  that  occur 
in  each  collision  [Refs.  1,2, 3,4,5]. 

The  interatomic  potential  function  determines  the  action 
of  the  collision  and  the  mechanics  of  the  collision  cascade. 
The  potentials  may  be  determined  theoretically  such  as  by 
performing  Hartree  Fock  calculations  or  semi-empirically , 
that  is,  derived  from  experimental  results.  This  method  can 
be  more  useful  because  the  theoretical  models  contain 
approximations  which  limit  their  accuracy.  Among  the 
potentials  developed  to  date  the  screened  Coulomb  potential 
has  been  particularly  successful.  If  there  were  no 
electrons  around  the  nuclei  of  atoms  the  potential  between 
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two  atoms  would  oe  a  simple  Coulomb  type.  The  presence  of 
electrons  changes  the  potential  making  it  more  complicated. 
The  screened  Coulomb  potential  accounts  for  the  electrons  by 
Introducing  an  additional  factor  called  the  screening 
function.  Many  variations  of  the  screening  function  have 
been  tried.  They  generally  produce  good  results  in  one 
energy  range,  but  poor  results  in  other  ranges.  More  than 
one  potential  can  be  pieced  together  with  splines  to  form  a 
continuous  potential  over  the  energy  range  of  interest. 
These  piecewise  potentials  usually  produce  good  results  in 
sputtering  calculations  [Ref.  6]. 

B.  LIQUIDS 

Much  sputtering  research  to  date  has  involved  the 
bombardment  of  solids  by  inert  gas  ions.  Solid  targets  can 
take  the  form  of  a  single  crystal,  a  polycrystal,  or  an 
amorphous  structure.  In  a  single  crystal  there  is  long 
range  order,  with  each  atom  located  at  a  lattice  point.  In 
an  amorphous  structure  the  atoms  still  maintain  short  range 
order  between  neighbors,  but  the  long  range  order  is 
eliminated. 

Liquids  are  similar  to  the  amorphous  solid,  in  that  they 
lack  long  range  order.  The  thermal  energy  associated  with 
the  atoms  of  the  liquid  allows  them  to  vibrate  further  from 
their  equilibrium  positions,  thus  permitting  the  possibility 
of  atoms  exchanging  positions  and  not  returning  to  their 
original  equilibrium  position.  This  weakening  of  the 
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crystal  structure  prevents  any  definition  of  a  fixed 
location,  and  therefore  the  material  has  no  rigid  structure 
[Ref.  7]. 

C.  HISTORICAL  OVERVIEW  OP  LIQUID  SPUTTERING 

Since  most  sputtering  research  has  been  done  with 
solids,  the  area  of  liquid  sputtering  is  in  a  relative 
infancy.  One  contribution  from  solid  target  research  which 
seems  applicable  to  liquids  is  that  for  low  energy  ion 
bombardment,  an  important  factor  effecting  the  yield  is  the 
surface  condition  of  the  target. 

In  1970,  Krutenat  and  Panzera  [Ref.  8]  looked  at  the 
yields  from  sputtering  liquid  and  solid  tin  with  Ar+  ions  of 
energies  from  threshold  to  1200  eV,  and  at  temperatures 
close  to  the  melting  point.  One  expected  difference  between 
the  two  surfaces  was  that  the  solid  would  build  up  the 
damage  from  the  ion  bombardment,  whereas  the  liquid  would 
"heal"  rapidly  from  previous  bombardments,  and  thus  present 
a  similar  target  to  all  incident  ions.  An  experimental 
measurement  of  the  yields  for  both  polycrystalline  solid  and 
liquid  tin  showed  that  for  ions  of  energy  up  to  about  400 
eV,  the  solid  had  a  higher  yield  than  the  liquid.  At  about 
400  ev  the  yields  were  the  same,  and  the  liquid  yield  was 
higher  above  400  ev.  Also,  the  yield  for  a  remelted  solid 
target  with  an  incident  ion  energy  of  about  475  ev  produced 
a  yield  that  was  between  the  liquid  and  solid  yields. 


Hurst  and  Cooper  [Ref.  9]  did  an  experiment  in  which 
they  measured  the  yield  from  an  indium  target  bombarded  with 
Ar+  ions  of  energy  107  eV,  as  they  raised  the  temperature  of 
the  target,  causing  it  to  melt.  This  gave  them  data  in  both 
the  solid  and  liquid  phases  for  the  same  ion  energy.  They 
also  ran  the  experiment  in  reverse  by  lowering  the 
temperature,  thus  changing  the  target  from  liquid  to  solid. 
They  found  that  the  liquid  produced  a  higher  yield  than  the 
solid  at  this  energy  which  is  close  to  threshold.  They  also 
ran  the  experiment  varying  the  ion  energy  from  17  ev  to  190 
eV  bombarding  both  solid  and  liquid  targets  at  temperatures 
close  to  the  melting  point,  and  found  that  the  liquid  yield 
exceeded  the  solid  yield  for  all  energies,  but  that  the 
yield  curves  had  the  same  shape  over  this  range.  They 
inferred  that  the  difference  in  the  yields  was  caused  by  the 
difference  in  binding  energies  of  the  two  surfaces,  but  the 
sputtering  effects  due  to  the  surface  condition  were 
approximately  the  same. 

In  1982  Dumke,  et  al.  [Ref.  10]  did  an  experiment  in 
which  they  sputtered  the  solid  and  liquid  phases  of  gallium, 
indium,  and  a  eutectic  gallium- indium  alloy,  using  argon 
ions.  They  found  no  significant  change  in  the  sputtering 
yields  between  the  solid  and  liquid  phases  of  the  pure 
elements.  The  solid  alloy  showed  a  surface  layer  of  indium 
which  dominated  the  sputtering  until  the  surface  layer  was 
sputtered  away  and  the  bulk  material  reached.  For  the 


liquid  alloy,  again  the  surface  was  mostly  indium,  even 
though  the  alloy  was  83.5%  gallium  and  16.5%  indium.  The 
indium  sputtering  yield  was  much  higher  than  for  gallium, 
indicating  that  most  sputtered  atoms  must  come  from  the 
surface  monolayer. 

D.  COMPUTER  SIMULATION 

Computer  simulation  of  physical  systems  is  now  being 
accepted  as  a  useful  adjunct  to  the  more  traditional 
theoretical  and  experimental  approaches.  Many  simulations 
of  sputtering  have  been  completed  over  the  thirty  years 
since  computers  first  became  available.  Few  simulations 
with  liquid  targets  have  been  done  to  date. 

D.  Y.  Lo,  et  al.  [Ref.  11]  investigated  collision 
cascades  in  liquid  Indium  using  a  multiple  interaction 
computer  simulation  code,  SPUT1.  The  liquid  indium  target 
was  melted  from  a  face  centered  cubic  (FCC)  structure  and 
heated  to  approximately  900  K.  The  calculations  were 
performed  using  pair  potentials.  A  Moliere  potential  was 
used  for  interactions  between  the  ion  and  target  atom.  For 
interactions  between  two  atoms  in  the  target  the  simulation 
used  a  potential  comprised  of  a  Moliere  core  joined  to  a 
Morse  well,  with  a  cubic  spline  [Ref.  11].  Analysis  of  the 
results  by  energy  and  angle  resolved  techniques  showed 
relatively  good  agreement  with  experimental  results  and  the 
theoretical  predictions  of  Garrison  [Ref.  11].  They 
concluded  that  the  angular  distribution  of  atoms  ejected 


from  the  first  layer  showed  qualitative  agreement  with 
experiment,  and  the  quantitative  disagreement  was  probably 
attributable  to  the  inadequacy  of  using  pair  potentials  in 
the  calculations. 
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II .  OBJECTIVES 


Most  of  the  sputtering  research  done  to  date  has  used 
the  method  of  pair  potentials  to  perform  the  calculations. 
First  a  potential  function  is  chosen  for  calculating  the 
potential  between  two  atoms  in  the  sample.  To  calculate  the 
force  on  a  given  atom,  all  the  atoms  in  the  sample  are 
paired  with  the  atom  located  at  the  point  in  question,  and 
the  potential  is  calculated  for  each  pairing.  The  summation 
of  the  potential  over  all  pairings  yields  the  total 
potential  of  the  atom  in  question.  By  repeating  this 
process,  the  potential  of  every  atom  in  the  sample  can  be 
determined.  From  the  potential  the  force  can  be  derived. 

The  pair  potential  method  has  limitations  on  its 
applicability  to  certain  physical  systems.  This  method 
gives  good  results  for  the  internal  collision  cascades,  but 
the  calculation  is  less  accurate  in  dealing  with  surface 
phenomena  [Ref.  12], 

A  new  method  of  calculating  the  potential  of  a  given 
atom  in  a  sample  was  proposed  in  1983.  This  new  method, 
called  the  embedded  atom  method  (EAM) ,  was  first  proposed  by 
Daw  and  Baskes  [Ref.  12].  In  the  embedded  atom  method  all 
atoms  in  the  sample  are  viewed  as  being  embedded  in  the  host 
consisting  of  all  other  atoms.  The  energy  depends  on  the 
electron  density.  This  overcomes  some  of  the  problems 
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associated  with  large  volume  dependent  terms  inherent  in  the 
pair  potential  approach.  The  semi-empirical  nature  of  the 
embedded  atom  approach  allows  the  potentials  to  be 
determined  in  terms  of  well  defined  and  measurable 
quantities . 

Both  methods  described  above  have  been  applied  to  liquid 
targets  to  a  limited  degree,  but  there  is  some  doubt  as  to 
the  validity  of  using  the  pair  potential  method  to  calculate 
the  energetics  of  an  amorphous  target.  On  the  other  hand 
the  EAM  avoids  some  of  the  drawbacks  of  the  pair  potential 
method,  and  can  still  handle  the  lack  of  structure 
associated  with  a  liquid  [Ref.  13].  For  this  reason  the  EAM 
will  be  used  for  all  the  calculations  in  this  project. 

The  main  objective  of  this  project  is  to  make  a  direct 
comparison  of  the  results  produced  by  calculating  the 
sputtering  yield  of  an  ion  bombardment  system  in  two 
different  states,  solid  and  liquid.  The  solid  target  will 
consist  of  a  perfect  crystal,  while  the  liquid  target  will 
be  the  same  crystal  which  has  been  warmed  above  its  melting 
point.  This  analysis  will  allow  conclusions  to  be  drawn  as 
to  the  effect  of  structure  or  long  range  order  on  a 
sputtering  yield. 

Different  types  of  trajectory  sets  for  the  liquid 
targets  will  be  used  in  order  to  investigate  the  effect  of 
sampling  methods  on  a  liquid  target  system. 
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All  of  the  calculations  and  results  derived  in  this 
project  will  be  accomplished  using  digital  simulation 
techniques.  In  order  to  ensure  a  proper  comparison,  all 
calculations  will  be  completed  using  the  same  target,  with 
the  exception  that  atoms  in  the  liquid  target  will  have  been 
displaced  from  their  crystal  lattice  positions,  an  amount 
which  on  average  is  equal  to  the  mean  thermal  displacement. 
Such  things  as  target  element,  size,  impact  points,  etc. 
will  remain  fixed. 


III.  MODEL  AND  SIMULATION 

A.  MULTIPLE  INTERACTION  MOLECULAR  DYNAMICS  SIMULATION 

The  computer  code  used  to  model  the  Ion  bombardment 
simulation  In  this  project  was  developed  by  Professor  Don  E. 
Harrison,  Jr.  at  the  Naval  Postgraduate  School.  The  program 
used  in  this  project  is  called  EWARM.  It  is  an  offshoot  of 
predecessors  called  EDYN  and  QDYN.  QDYN  is  the  grandfather 
of  all  of  them.  It  had  been  developed  to  do  sputtering 
simulation  research  using  pair  potentials.  EDYN  is  similar 
to  QDYN  except  that  it  uses  the  EAM  to  calculate  the  forces 
needed  for  the  simulation.  A  detailed  account  of  the 
digital  simulation  method,  the  mathematical  model  and  the 
operation  of  the  integration  technique  has  been  previously 
published  by  Harrison,  et  al.  [Refs.  14,  15,  16,  17,  18, 

19]- 

EWARM  differs  from  the  other  two  in  that  it  contains  the 
ability  to  handle  the  simulation  of  a  liquid  target.  The 
major  difference  is  in  the  addition  of  a  WARMER  subroutine 
which  randomly  dislocates  the  atoms  of  the  target  from  their 
perfect  crystal  lattice  sites,  thus  producing  an  amorphous 
target . 

All  three  of  these  programs  handle  the  molecular 
dynamics  in  the  same  way.  EDYN  is  completely  contained 
within  EWARM.  EWARM  is  a  continuous  time  simulation  of 
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sputtering,  which  means  that  the  program  proceeds  through  a 
series  of  timesteps.  Each  timestep  consists  of  a 
calculation  of  the  forces,  the  new  velocities  and  the 
positions  of  each  atom,  at  the  end  of  the  timestep,  movement 
of  the  atoms  to  their  new  positions,  and  a  test  for  energy 
conservation.  Both  programs  compute  the  atoms'  trajectories 
by  performing  a  numerical  solution  of  Newton's  classical 
equations  of  motion.  The  length  of  subsequent  timesteps  is 
controlled  by  a  predictor-corrector  integration  scheme  which 
determines  the  time  increment  from  the  fastest  moving 
particle  [Ref.  17]. 

The  system  used  for  this  project  consists  of  an 
impinging  argon  ion  and  a  target  of  912  rhodium  atoms.  Each 
particle  in  the  system  is  characterized  by  mass,  position, 
velocity,  and  the  force  laws  by  which  it  interacts  with 
other  particles. 

After  the  impinging  ion  hits  the  target  at  the  specified 
impact  point,  the  trajectories  of  all  the  particles  in  the 
system  develop  in  time  as  the  energy  and  momentum  are 
dissipated  through  the  target.  The  calculation  is  called  a 
trajectory  and  the  dissipation  through  the  target  is  called 
a  collision  cascade.  The  trajectory  is  terminated  when  the 
energy  and  momentum  has  dissipated  to  the  point  where  no 
more  atoms  will  be  ejected.  in  this  particular  system  the 
calculation  ends  when  the  maximum  kinetic  energy  of  a  single 
particle  falls  below  2.5  ev. 
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B.  THE  EMBEDDED- ATOM  METHOD 


In  the  embedded -atom  method  each  atom  is  treated  as  if 
it  were  embedded  in  the  host  of  all  other  atoms.  The  energy 
required  to  embed  an  atom  depends  on  the  electron  density. 
This  electron  density  can  always  be  defined ,  so  any  problems 
associated  with  not  being  able  to  define  a  volume  are 
avoided.  Therefore  the  embedded  atom  method  can  be  used  on 
targets  with  surfaces,  and  promises  better  results  when 
applied  to  liquid  targets.  Another  important  factor  is  that 
the  embedded-atom  method  in  not  significantly  more 
complicated  to  use  than  pair  potentials  [Ref.  12]. 

To  begin  we  think  of  each  atom  as  an  impurity  embedded 
in  the  host  material.  The  total  potential  energy  is  the  sum 
of  the  host  and  impurity  contributions.  The  host  potential 
is  a  function  of  the  electron  density  without  the  impurity, 
and  the  impurity  potential  depends  on  the  position  and 
charge  of  the  impurity  nucleus.  Therefore  the  energy  of  the 
host  with  impurity  is  a  functional  of  the  host  electron 
density  and  a  function  of  the  impurity  type  and  position. 
The  functional  is  universal,  independent  of  the  host 
material.  Its  form  is  not  known,  but  a  simple  approximation 
assumes  that  energy  only  depends  on  the  other  atoms  that  are 
close  to  the  impurity,  or  the  impurity  experiences  a  locally 
uniform  electron  density.  with  this  approximation  the 
functional  of  the  electron  density  becomes  a  function  and 
only  the  electron  density  at  the  impurity  position  is 
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needed.  Th* *?  makes  the  position  dependence  trivial  and  the 
energy  depends  only  on  the  electron  density  of  the  host  at  a 
given  point  plus  an  electrostatic  interaction.  The  total 
energy  is  a  sum  over  all  contributions  from  the  individual 
atoms.  The  electrostatic  interaction  and  the  electron 
density  terms  are  determined  empirically  [Refs.  10,  12]. 
Daw  and  Baskes  have  previously  published  articles  describing 
the  approximations  and  details  of  the  derivation  of  the 
embedded  atom  method  [Refs.  10,  12]. 

A  new  version  of  QDYN  was  created  using  the  embedded- 
atom  method  to  calculate  the  potentials  of  the  atoms  in  the 
system.  Called  EDYN,  it  employs  a  method  whereby  the  total 
potential  energy  of  a  given  atom  is  the  sum  of  the  short- 
range  pairwise  potential  energy  and  the  embedding  energy. 

C.  WARMING  THE  TARGET 

The  next  modification  implemented  in  QDYN  and  EDYN  was 
the  addition  of  the  subroutine  WARMER.  This  routine  takes  a 
perfect  crystal  lattice  and  randomly  moves  the  atoms  so  that 
the  resulting  target  is  amorphous.  Specifically,  this  is 
done  by  adding  a  correction  term  to  each  component  of  the 
initial  positions  of  the  perfect  crystal  target  atoms.  This 
correction  term  is  simply  a  normally  distributed  random 
variant  multiplied  by  a  thermal  amplitude  factor.  The 
temperature  can  be  approximated  from  the  change  in  the 
potential  energy  of  the  system  so  that  the  target  can  be 
"warmed”  to  a  desired  temperature.  If  the  target  is  wanned 
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to  the  melting  temperature  of  the  material,  then  the 
resulting  target  is  a  liquid. 

Since  the  warmer  randomly  locates  the  atoms  of  the 
target,  it  needs  a  supply  of  random  numbers.  These  come 
from  a  subroutine  that  generates  pseudo-random  numbers.  As 
with  most  pseudo-random  number  generators,  the  user  has 
control  of  the  starting  point  of  the  generator  through  the 
seed.  Each  number  produced  by  the  generator  is  then  used  as 
the  seed  for  the  next  number.  Throughout  this  project  the 


starting  and  ending  seeds  are  monitored  so  that  a 
calculation  can  be  repeated  by  using  the  same  starting  seed. 

These  starting  seeds  are  very  important  because  changing 
the  seed  will  change  the  way  in  which  the  target  atoms  are 
dislocated,  and  therefore  produce  a  different  target.  For 
example,  to  run  a  set  of  trajectories  with  the  same  target 
would  require  that  the  same  starting  seed  be  used  for  the 
warming  of  each  trajectory  target.  This  line  of  thinking 
can  lead  to  many  possible  ways  of  warming  a  target  for  a  set 
of  trajectories. 

D.  POTENTIAL  FUNCTIONS 

In  a  simulation  of  this  type,  the  interatomic  potential 
functions  must  be  chosen  so  that  the  forces  between  atoms 
can  be  calculated.  These  functions  ultimately  control  the 
dynamics  of  the  simulation.  Approximations  will  occur, 
because  of  a  lack  of  detailed  knowledge  in  certain  regions 
of  the  interatomic  potentials.  At  separations  much  greater 


wr  ^ 


than  the  equilibrium  separation  of  two  atoms,  there  is 
relatively  good  understanding  of  the  potentials.  Likewise, 
at  the  equilibrium  separation,  but  at  the  separation 
distances  less  than  equilibrium,  and  in  the  transition  area 
where  the  separation  distance  is  just  greater  than 
equilibrium,  much  less  is  known. 

There  are  several  theories  available  as  well  as  some 
empirical  methods  for  obtaining  potential  functions  to  use 
in  a  calculation.  Most  of  these  potential  functions  take 
the  form  of  a  "wall"  or  almost  vertical  slope  as  the 
interatomic  distance  approaches  zero,  decreasing  to  a  "well" 
or  minimum  at  the  equilibrium  separation,  increasing  through 
a  steady  transition  approaching  a  horizontal  line  far  from 
the  equilibrium  separation  distance.  This  form  of  a 
potential  function  will  produce  an  attractive  force  at 
distances  greater  than  the  equilibrium  separation,  zero 
force  at  equilibrium,  and  a  repulsive  force  which  gets  very 
large  in  the  area  where  the  separation  distance  approaches 
zero  [Ref.  6]. 

In  most  cases,  one  function  does  not  describe  the 
interatomic  potential  very  well  over  all  of  the  regions 
described  above.  Several  functions  are  generally  splined 
together  to  produce  a  composite  function  which  is  good  over 
all  regions,  and  can  be  verified  by  existing  experimental 
data. 
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In  this  project  the  interatomic  potentials  are 
empirical.  They  were  created  from  several  functional  forms 
which  agree  with  experimental  data  in  the  region  where  they 
are  employed.  They  are  functions  that  use  certain  constants 
related  to  physical  parameters  of  the  atoms  involved. 

The  functions  used  in  this  project  to  describe  the 
interatomic  potentials  can  be  divided  into  two  types.  The 
first  describes  the  interactions  between  the  ion  and  a 
target  atom  (Ar-Rh),  and  the  second  describes  the 
interactions  between  two  target  atoms  (Rh-Rh) . 

The  Rh-Rh  function  consists  of  a  modified  Moliere 
potential  joined  to  an  attractive  Morse  potential  using  a 
cubic  spline.  This  is  given  by: 


v  -  [ (Z1Z2e2/ka)/(R/ka) ]  g(R/ka) 


R  <  R, 


g(R/ka)  -  [0.35  exp(-0.3  R/ka)  +  0.55  exp(-1.2  R/ka) 

+  0.10  exp(-6.0  R/ka)] 


a  -  O.easSao/tZi1/2  +  Z21/2]2/3 
V  -  c0  +  CXR  +  C2R2  +  C3R3 


Ra  <  R  <  Rb 


V  -  D{exp[-2b(R-Re) ]  -  2  exp[-b(R-Re) ] } 


Rb  <  R  <  Rc 


V  -  0 


R  >  Rr 


where : 


V  is  the  potential 

Zi  is  the  atomic  number  of  the  first  atom 
Z2  is  the  atomic  number  of  the  second  atom 
e  is  the  charge  of  an  electron 
R  is  the  distance  between  the  two  atoms 
k  is  a  magnitude  factor  between  0  and  l 


a.l  . 
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ag  is  the  Bohr  radius  (0.529  A) 
g  is  the  screening  function 
Cg,  Ci,  C2,  C3  are  fitting  parameters 
D  is  the  well  depth 

Re  is  the  equilibrium  separation  distance 
b  is  a  fitting  parameter  (A-1) 


In  reality  the  potential  function  exists  out  to  infinite 
separation.  The  composite  function  is  truncated  at  Rc 
because  beyond  this  point  its  effect  is  negligible  [Ref.  6]. 

The  potential  function  that  was  used  for  interactions 
between  Ar+  and  Rh  is  an  unmodified  Moliere  function  of  the 
form: 


V  -  [ (Z1Z2e2/a)/(R/a) ]  g(R/a) 

V  -  0 


R  <  R, 


R  >  R, 


The  parameters  that  define  the  potential  functions  for 
all  the  runs  generated  in  this  simulation  are  given  in 
table  l. 


TABLE  1 

POTENTIAL  FUNCTION  PARAMETERS 


Parameter 

Rh-Rh 

Ar-Rh 

D  (ev) 

0.7595 

0.0 

Re  (A) 

2.750 

0.0 

b  (A'1) 

1.080 

0.0 

Ra  (LU) 

0.73 

1.70 

Rb  (LU) 

0.80 

1.70 

Rc  (LU) 

2.05 

1.70 

k 

0.74 

1.0 

E .  PHYSICAL  PARAMETERS 

The  target  is  produced  by  building  a  perfect  crystal, 
and  then  "warming"  the  crystal  to  the  melting  point  of 
rhodium  (2239  K),  thus  producing  an  amorphous  liquid.  The 
target  is  oriented  so  that  the  (111)  face  of  the  perfect 
crystalline  target  would  be  normal  to  the  incoming  ion. 
This  notation  of  a  (111)  face  is  not  relevant  for  a  liquid 
but  does  give  a  basis  for  comparison  with  results  from  the 
solid  crystal. 

Rhodium  is  a  face  centered  cubic  metal  (fee)  with  a 
lattice  constant  (a0)  of  3.804  A.  The  basic  unit  of 
distance  used  is  the  lattice  unit  which  is  defined  as  half 
the  interatomic  spacing.  For  this  simulation  the  lattice 
unit  is  equal  to  1.902  A. 

In  the  simulations  it  is  important  to  choose  the  target 
size  judiciously.  A  target  size  must  be  chosen  that  will 
produce  realistic  results  yet  keep  within  the  constraints  of 
the  computer  resources  available. 

After  extensive  experience  with  these  types  of 
simulations,  Harrison  has  developed  a  working  definition  of 
containment  [Ref.  17]:  "If  an  increase  in  target  size  does 
not  change  the  results  of  the  computation,  the  trajectory  is 
effectively  contained  for  the  purposes  of  that  computation." 

As  the  target  size  increases,  the  observable  will 
approach  its  absolute  value  asymptotically  as  full 
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containment  is  reached.  In  further  clarification  of  this 
point  Harrison  states  that: 

As  one  might  anticipate  from  experimental  experience 
with  absolute  values,  yield  provides  the  most  severe 
test  of  containment;  so  absolute  yield  containment 
guarantees  everything  else.  As  the  target  size 
increases,  all  other  global  results  reach  constant 
values  for  much  smaller  targets  than  those  required  to 
produce  stable  yield  values  [Ref.  17]. 

F.  ENSEMBLES  OF  TRAJECTORIES 

Another  item  requiring  careful  choice  is  the  method  by 

which  the  trajectories  are  chosen  so  that  the  result  is  an 

average  of  the  observable  that  we  wish  to  measure.  A  method 

has  been  developed  for  perfect  crystal  targets,  and  has  also 

been  applied  to  warmed  targets  [Ref.  19].  The  method  uses 

the  concept  of  a  representative  area.  Each  impact  point  on 

the  target  surface  produces  a  different  trajectory.  By 

symmetry,  a  subset  of  impact  points  on  the  target  surface 

will  produce  every  possible  trajectory  for  a  given  lattice 

orientation.  The  representative  area  is  this  subset  of 

impact  points  that  contains  every  possible  trajectory  for  a 

given  target  surface.  Since  there  are  an  infinite  number  of 

possible  trajectories,  the  sample  size  requires  judicious 

choice.  A  carefully  planned  sample  of  trajectories  over 

this  representative  area  should  produce  the  average  value 

for  the  observable.  The  set  of  all  possible  trajectories  is 

called  an  ensemble,  and  the  smaller  sample  set  is  known  as  a 

sample  ensemble.  Regular  sampling  of  the  representative 

area  using  upwards  of  300  impact  points  has  been  used  very 


successfully  in  the  past  and  will  be  used  in  this 
simulation.  A  possible  representative  area  is  a  rectangle 
for  the  (ill)  surface.  A  starting  point,  called  a  pinpoint 
is  placed  in  the  lower  left  corner  of  the  representative 
area,  and  a  set  of  impact  points  are  regularly  spaced  from 
the  pinpoint  over  the  representative  area  [Ref.  19]. 

For  this  project  a  set  of  300  impact  points  has  been 
chosen.  This  set  of  impact  points  will  be  used  for  the 
perfect  crystal  target  as  well  as  for  the  liquid  targets. 
The  use  of  this  sampling  method  has  been  substantiated 
previously  for  perfect  crystal  targets  [Ref.  19],  but  when 
the  target  melts  the  justification  for  its  use  is  more 
questionable.  This  project  aims  to  determine  a  legitimate 
method  of  sampling  for  liquid  targets. 

G.  ANALYSIS 

The  data  from  the  simulation  program  can  be  analyzed 
using  a  fortran  program  called  ANPLOT  which  takes  the 
information  about  the  individual  ejected  atoms  and  produces 
graphical  presentations  allowing  comparison  to  other 
simulation  or  experimental  results.  It  produces  various 
plots  and  histograms  such  as  atoms  per  single  ion,  ejected 
atom  energy  distribution,  atom  ejection  time  distribution, 
etc.  This  program  was  used  to  analyze  all  the  data  derived 
from  this  project. 
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IV.  RESULTS 

A.  OVERVIEW 

The  results  presented  in  this  project  were  produced  by  a 
series  of  simulations.  The  total  amount  of  computer  time 
used  to  collect  this  information  is  in  excess  of  90  hours. 
A  single  sampling  method  was  employed  for  the  solid  case, 
but  six  different  methods  were  tried  with  liquid  targets  in 
order  to  investigate  sampling  methods  to  be  used.  Of  the 
six  methods  tried,  four  appear  to  be  useful  for  sputtering. 
In  the  solid  case  300  trajectories  were  run  with  a  perfect 
crystal  target.  The  target  was  the  same  for  each 
trajectory.  The  first  liquid  set  consisted  of  300 
trajectories  hitting  a  target  that  had  been  wanned.  The 
target  was  identical  for  each  trajectory.  The  next  data  set 
also  consisted  of  300  trajectories,  but  a  different  seed  was 
used  to  produce  the  target  for  each  trajectory,  so  that  each 
trajectory  was  hitting  a  different  target.  Two  sets  of  100 
trajectories  were  run  in  which  the  impact  points  remained 
fixed  while  the  seed  was  changed,  thus  producing  different 
targets  for  each  trajectory. 

For  the  last  three  sets,  four  atoms  in  the  initial 
impact  area  were  frozen  in  their  crystal  lattice  positions 
while  all  the  other  atoms  were  randomly  dislocated  by  the 
warmer.  These  four  atoms  enclose  the  representative  area  in 
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the  perfect  crystal  target,  and  so  it  was  postulated  that  by 
freezing  them,  some  credence  could  be  given  to  using  this 


sampling  technique  on  a  liquid  target.  One  set  of  300 
trajectories  was  run  with  the  same  warming  seed  for  each 
trajectory,  one  set  of  300  trajectories  run  with  a  different 
seed  for  each  trajectory,  and  finally  the  last  set  of  100 
trajectories  was  run  keeping  the  impact  point  fixed  and 
changing  the  seed  for  each  trajectory. 

The  trajectory  sets  completed  with  fixed  impact  points 
were  analyzed  as  distributions  to  determine  the  possibility 
of  using  this  method  as  a  sampling  technique,  and  to  see 
what  kind  of  variation  would  be  produced  at  a  single  impact 
point.  From  all  the  data  it  is  apparent  that  changing  the 
seed  does  effect  the  trajectory,  and  therefore  the 
sputtering  yield. 

B.  SPUTTERING  YIELD 

The  total  sputtering  yields  calculated  for  the  four  300 
trajectory  sets  using  liquid  targets  were  approximately  the 
same.  There  is  less  than  a  2%  difference  between  any  single 
set  average  and  the  average  of  the  four.  The  fact  that 
these  four  sets  independently  reproduced  the  same  relative 
yield  would  indicate  that  the  sample  is  probably  large 
enough  for  this  observable.  It  does  not  give  any  indication 
that  any  of  the  sampling  methods  is  better  than  another,  and 
this  paper  has  not  tried  to  optimize  the  sample  size  or 
method.  The  average  of  these  four  sets  indicates  that  the 
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liquid  target  produced  about  a  6%  higher  yield  than  the 
solid  target.  The  single  impact  point  trajectory  sets 
showed  reproducibility,  but  did  not  agree  with  the  other 
liquid  sets.  This  means  that  for  our  method  of  warming  the 
crystal,  changing  the  impact  point  on  the  target  is  not  the 
same  as  changing  the  target  behind  the  impact  point, 
probably  due  to  the  relatively  small  thermal  displacements 
of  the  atoms  from  their  initial  positions.  This  should  be 
explored  further  in  any  future  work.  The  total  sputtering 
yields  are  presented  in  table  2. 


Category 


TABLE  2 

TOTAL  SPUTTERING  YIELDS 

Atoms  Average  Atoms 

Ejected  Trajectories  per  Trajectory 


perfect  crystal 

1596 

300 

5.32 

liquids: 

single  target 

1670 

300 

5.57 

different  targets 

1731 

300 

5.77 

liquids  with  4  representative  area  atoms 

initially  f: 

single  target 

1706 

300 

5.69 

different  targets 

1725 

300 

5.75 

liquids  with  fixed 
point  l 

impact  point 
795 

102 

7.79 

point  2 

748 

103 

7.26 

point  1* 

803 

108 

7.44 

same  impact  point,  except  that  four  atoms  in  the 
representative  area  are  fixed 
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Because  of  the  apparent  difference  between  the  single 
impact  point  sets  and  the  other  liquid  sets,  the  multiple 
impact  point  sets  will  be  compared  to  the  solid  set  in  the 
remainder  of  the  results.  The  single  impact  point  sets  are 
too  locally  controlled  by  the  choice  of  impact  point.  The 
multiple  impact  point  sets  give  a  better  sampling  diversity. 

C.  ATOMS  EJECTED  PER  SINGLE  ION 

Atoms  per  single  ion  (ASI)  is  an  important  quantity 
which,  in  sputtering  simulations  of  solids  can  be  determined 
relatively  accurately  using  the  molecular  dynamics 
simulation  method.  The  sputtering  yield  is  the  mean  ASI 
over  the  sample  set.  The  ASI  distribution  is  difficult  to 
obtain  in  an  experiment,  and  simulations  can  often  gain 
useful  physical  information.  In  this  simulation  the  liquid 
ASI  distributions  did  not  differ  significantly  from  the 
solid  ASI  distribution,  except  in  producing  a  higher  mean 
(Table  2).  Figures  1  through  5  show  the  ASI  distributions 
for  the  solid  trajectory  set  and  the  four  liquid  trajectory 
sets.  The  two  liquid  distributions  (Figures  2,3)  that  used 
the  same  target  for  all  trajectories  closely  follow  the 
solid  distribution  (Figure  1).  The  other  two  liquid  targets 
that  used  a  different  seed  for  each  trajectory  (Figures  4,5) 
only  differ  slightly. 

D.  EJECTION  TIMES 

Figures  6  through  10  provide  a  view  of  the  ejection  time 
distributions  for  the  five  cases  of  interest.  There  is  no 
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apparent  difference  in  the  ejection  time  distributions 
between  the  liquid  and  solid  cases.  Also,  the  sampling 
method  appears  to  make  little  difference  in  the  liquid 
cases.  An  important  point  is  that  about  95%  of  all  the 
ejections  of  this  system  occur  within  first  200  femtoseconds 
of  the  simulation,  much  shorter  than  the  time  scale  for 
thermal  effects. 

E.  EJECTED  ATOM  ENERGY  DISTRIBUTIONS 

There  is  very  little  difference  between  any  of  the 
energy  distributions  of  the  ejected  atoms,  meaning  that  the 
state  of  the  target  does  not  effect  the  energy  distribution 
of  the  sputtered  particles.  Also,  the  sampling  method  does 
not  appear  to  effect  the  resulting  energy  distributions. 

F.  EJECTED  ATOM  ANGULAR  DISTRIBUTIONS 

Sputtering  simulations  can  also  determine  the  angular 
distributions  of  the  atoms  as  they  are  ejected.  This  is  a 
simulation  observable  that  compares  well  with  experimental 
sputtering  data  [Ref.  20].  A  set  of  plots  has  been 
developed  by  Harrison  to  show  the  polar  angular 
distributions.  These  distributions  are  created  by  tallying 
the  ejected  atoms  into  bins  five  degrees  wide  from  0  to  90 
degrees.  0  degrees  is  normal  to  the  surface.  If  the  low 
energy  atoms ,  say  below  five  ev  are  excluded,  the  comparison 
to  the  experimental  data  is  even  better.  In  this  project 
all  of  the  ejected  atoms  were  included  in  the  distributions. 
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The  distributions,  presented  as  figures  ll  through  15, 
show  little  difference  between  the  liquid  and  solid  cases. 
Also,  there  is  no  detectable  difference  among  the  four  cases 
of  the  liquid.  The  peaks  of  these  distributions  all  occur 
about  30  degrees. 

By  excluding  those  ejected  atoms  with  less  than  20  ev  of 
energy,  the  distributions  are  skewed  to  smaller  polar 
angles,  meaning  that  the  higher  energy  atoms  must  be  ejected 
at  angles  closer  to  normal  than  those  atoms  with  lower 
energy.  Figure  16  is  provided  as  a  comparison  for  the  solid 
case.  The  liquid  distributions  show  the  same  trend. 

G.  SPOT  PATTERNS 

Spot  patterns  are  another  method  of  showing  the  angular 
distributions  of  the  ejected  atoms.  Spot  patterns  have  an 
advantage  over  the  polar  angular  distribution  plots,  namely 
that  the  spot  patterns  present  the  azimuthal  angular 
dependence  simultaneously.  Figures  17  through  19  show  the 
spot  patterns  for  the  solid  case  and  the  four  liquid  cases. 
The  solid  case  shows  its  distinctive  six  fold  symmetry 
pattern.  In  the  liquid  cases  this  pattern  is  still  evident, 
but  not  as  sharp. 

H.  ATOM  YIELD  PER  IMPACT  POINT 

Figures  20  through  24  are  the  graphical  representations 
of  the  relative  yields  produced  at  each  impact  point  for  the 
five  cases  being  studied.  This  displays  each  impact  point 
on  the  target  surface,  and  then  indicates  the  relative  yield 
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at  that  point  by  a  series  of  tree  rings.  Higher  yields  are 
shown  by  more  rings.  The  locations  of  the  centers  of  the 
first  surface  layer  target  atoms  have  been  plotted  to  give  a 
better  picture  of  where  the  atoms  are  relative  to  the  impact 
points.  In  the  three  targets  in  which  the  seed  was 
identical  for  each  trajectory  (Figures  20,21,22),  the 
ejections  were  less  uniformly  distributed  than  for  those  in 
which  the  target  was  changed  from  trajectory  to  trajectory 
(Figures  23,24).  In  neither  case  were  the  ejections 
uniformly  distributed  over  the  representative  area,  which 
should  be  the  case  for  a  purely  random  target. 
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V.  CONCLUSIONS 

This  project  has  made  no  attempt  to  compare  the  results 
of  the  simulations  with  experimental  results.  The 
comparisons  reported,  have  been  between  similar  systems  in 
two  different  states,  solid  and  liquid,  and  between 
identical  models  of  liquid  systems  in  which  the  sampling 
methods  were  varied. 

All  of  the  sampling  techniques  used  produced  about  the 
same  results  for  the  calculated  yields.  Changing  the  seed 
to  the  random  number  generator  changes  the  trajectory,  as 
shown  by  the  distributions  produced  at  the  same  spatial 
impact  point.  Figure  25  is  an  example  of  the  distributions 
that  resulted  from  this  sampling  technique.  The 
distributions  for  the  other  observables  were  similar. 
Changing  the  impact  point  produced  better  statistics  than 
keeping  the  impact  point  fixed  and  changing  the  seed.  This 
is  because  the  thermal  displacements  were  with  respect  to 
positions  based  on  a  crystal  lattice.  Keeping  the  impact 
point  fixed  maintains  some  of  the  underlying  crystal 
structure,  and  hence  the  target  is  not  truly  randomized. 
Computer  simulation  of  sputtering  from  liquids  should  not 
use  this  sampling  technique. 

Any  of  the  other  four  methods  should  give  reasonable 
statistics.  It  appears  that  the  size  of  the  representative 
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area  chosen  for  solid  simulations  because  of  lattice 
symmetry  also  works  well  for  the  randomized  target. 

Another  alternative  is  to  use  a  stochastic  sampling 
technique  of  the  representative  area,  in  addition  to 
randomizing  the  target  for  each  trajectory. 

A  possible  improvement  to  the  model  would  be  to  simulate 
a  warmed  target  by  adding  thermal  velocities  in  addition  to 
displacements.  Any  initial  kinetic  energy  an  atom  may  have 
could  affect  its  ejection  probability.  It  is  recommended 
that  this  approach  be  implemented  in  any  future  studies  of 
liquid  targets. 

The  effect  of  warming  a  target,  even  past  its  melting 
point,  does  not  substantially  change  the  calculated 
sputtering  yields.  A  slight  increase  in  yield  occurs  as  a 
result  of  the  warming,  but  the  ejected  particle 
distributions  do  not  change. 

The  work  described  in  this  project  used  the  embedded 
atom  method  to  perform  all  the  calculations.  This  is  the 
new  and  promising  method,  but  most  current  research  still 
uses  pair  potential  calculations.  Other  many  body 
potentials  have  now  been  reported  for  modelling  equilibrium 
phenomena  [Ref.  21]  which  promise  to  give  even  better  models 
of  the  underlying  physics.  An  extension  of  this  project 
would  be  to  run  the  simulations  using  a  variety  of 
potentials,  and  to  compare  with  results  obtained  by 
experiment. 
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Further  study  of  different  sampling  techniques  would 
only  improve  the  confidence  in  any  sampling  method 
|  ultimately  selected . 
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APPENDIX 
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Figure  1  Atoms  per  Single  Ion  Distribution. 
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Figure  3  Atoms  per  Single  Ion  Distribution. 
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Figure  4  Atoms  per  Single  Ion  Distribution. 
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Figure  6  Atom  Ejection  Time  Distribution. 
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Figure  7  Atom  Election  Time  Distribution 
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Figure  8  Atom  Ejection  Time  Distribution. 


40 


1.0  RH(111)/AR<U1>  WRMD  DIFF  TARG  W/4  FROZEN 

ATOM  EJECTION  TIME  DISTRIBUTION 


500.0 


EJECTION  TIME  (FS) 


Figure  10  Atom  Ejection  Time  Distribution 
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Figure  13  Ejected  Atom  Angular  Distribution  -  All  Atoms 
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Figure  14  Ejected  Atom  Angular  Distribution  -  All  Atoms 
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Figure  15  Ejected  Atom  Angular  Distribution  -  All  Atoms 
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Figure  23  Atom  Yield  per  impact  Point. 
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Figure  25  Atoms  per  Single  Ion  Distribution. 
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